Mechanistic model of the Oxygen isotope in tree ring cellulose

Isotope Dendrochronology

COURSE: DISC, LTRR, UNIVERSITY OF ARIZONA. INSTRUCTORS: P. SZEJNER & S. BELMECHERI

Instructions:

The following code and comments should further your experience of using a mechanistic model to understand the isotopic fractionation mechanisms underlying Oxygen isotope (d18O) in tree ring cellulose. This exercise will focus on the main determinants (variables and parameters) of cellulose d18O including: (1) source water isotope ratio, (2) leaf water isotopic enrichment, and (3) biochemical fractionation at the site of sucrose production and cellulose synthesis.

The exercise and the mechanistic model correspond to the process based models termed “Craig-Gordon” model [1], and the “Roden” model [2].

[1] Craig H, Gordon L (1965) Deuterium and oxygen 18 variations in the ocean and the marine atmosphere. In: Tongiorgi E (ed) Stable Isotopes in Oceanographic Studies and Paleotemperatures. Consiglio Nazionale delle ricerche Laboratorio di Geologia Nucleare. Pisa, Spoleto, Italy, pp 9–130).

[2] Roden J, Lin G, Ehleringer JR (2000) A mechanistic model for interpretation of hydrogen and oxygen isotope ratios in tree-ring cellulose. Geochim Cosmochim Acta 64:21–35. https://doi.org/10.1016/S0016-7037(99)00195-7

The stable Oxygen isotope composition of cellulose (δ18Ocell) is determined by three factors/processes:
1. The isotopic content of the source water.
2. The leaf evaporative enrichment effect.
3. The proportion of phloem sugars that exchange with source water during cellulose formation.

The signal of these three processes is integrated into the cambium during xylogenesis. The goal of this exercise is to: (1) implement a mechanistic model to account for all these factors processes and their associated fractionations, (2) to use this mechanistic model to interpret variations of cellulose d18O.

You will simulate three different isotopic ratios in tree ring cellulose for different environmental conditions. These environmental condition reflect the growing season variability of relative humidity and precipitation experienced by trees growing in the US Southwest (See Figure 1). The growing season in this region is characterized by a bimodal distribution of precipitation (Winter-snowpack and Summer-Monsoon). The relative humidity varies from dry spring to humid summer. The two periods are separated by a hyper arid period where with low precipitation and extremely low relative humidity. Accordingly, the variables in the model will vary as follows:
* The relative humidity from 50% (spring) to 30% (pre-summer drought) to 60% (monsoon).
* The mean air temperature will vary from 25, to 33, to 27 Celsius.
* The source water will remain constant at -8 permill.

At the end of the exercise, you will compare the simulated 18O cellulose for the 3 seasonal variations of RH (and temperature) with measurements of cellulose d18O in tree rings from a site located in the US Southwest.

Figure 1: Seasonal relative humidity in HEL-Catalina Mountains (Arizona) (Belmecheri et al. 2018 PCE). Notice the large variations of relative humidity.

Instructions

You will need to download the R project from https://github.com/paulszejner/DISCO-Notebook2.git and set your own working directory. Once you open the Notebook2.Rproj in R studio you can open the Notebook2.Rmd file to proceed with this exercise.

PART 1: Source Water

The d18O of source water noted as Source_18O is the isotopic ratio of the water transported through the xylem up to the leaves. This is derived from direct precipitation, soil water or Thus, it is the starting point of d18O that will be modified by the isotopic effects (fractionations) during all the steps leading to its fixation in tree ring cellulose. Having a robust estimation or measurement of source water d18O is therefore paramount.

Here, the source water is set up at a fixed value:

Source_18O <- -8 # This is the isotopic value of the source water and  will remain constant at -8 permill during this  exercise

Then, it is converted to an isotopic ratio with regard to the standard VSMOW:

rsmow <- 0.0020052 # Oxygen Isotope ratio from STANDARD VSMOW
Rsource= rsmow*(1+Source_18O/1000) 

Next, the isotopic value of the water vapor is estimated based on source water d18O and leaf temperature, assuming equilibrium between atmospheric vapor and source water. Water vapor can be estimates as follows:

water_vapor = Source water + the equilibrium fractionation between liquid water and vapor.

The equilibrium fractionation varies as a function of temperature.

# water_vapor = Source water + the equilibrium fractionation.

# The 3 temperature values need to be modified to reflect the seasonality described in the introduction
temperature <- c(40, 40, 40) # temperature here is in Celsius 

#**** WARNING**** for water vapor temperature is in Kelvin (273.15)
water_vapor <-Source_18O +(-1*(1.137*(1000000/((temperature+273.15)^2))-0.4156*((1000/(temperature+273.15))-2.0667))) 

Then, it is converted to an isotopic ratio with regard to the standard VSMOW:

Rwater_vapor= rsmow*(1+water_vapor/1000)
#note!  I need help on figuring out this! What is needed?

PART 2: leaf water isotopic enrichment

During plant transpiration, the heavier H218O evaporates and diffuses more slowly through the stomata than does H216O, leaving leaf water enriched in 18O. At steady state and considering the leaf as a single and isotopically well-mixed water pool, the leaf water isotopic enrichment can be estimated using the Craig Gordon model.

#Constants

# br = boundary layer conductance. (mol m-2 s-1)
br <-1 #considered reasonable value between 0.02 to 3 in Roden paper boundary layer resistance    
ak <- 1.0285    # alpha_kappa, Kinetic Fractionation 
akb <- 1.0189   # alpha_kappa_beta, Kinetic Fractionation boundary layer


# Variables  remember to modify these values to the variability form the Figure 1

###########
Relative_humidity <- c(100, 100, 100)
###########

Stomatal_conductance <- 0.5

We already know that the source water has this values:

round (Rsource, digits = 4)

We already know that the water vapor could be close to has this values:

round(Rwater_vapor, digits = 4)

# a*  **equilibrium fractionation at leaf temperature  
# The equilibrium fractionation factor, a* is temperature dependent 
# (see Majoube, 1971, J. Chim. Phys. 58:1423-1436)
alpha_star =exp((1.137*1000/(273.15+temperature)^2)-(0.4156/(273.15+temperature))-0.0020667)
#R source O  isotope ratio, source water


# barometric    pressure  (mmbar)  
press_bar <- 77 #a good site average is 77.5 according to W, This need to be changed for the pressure level

Here we assume that the vapor is saturated inside the leaf. We apply this formula to calculate the vapor pressure when saturated (Rh=100%) at the temperature we best estimate the plant is and call it "ei".

ei <- (101325*exp((((-0.1299*(1-(373.16/(273.15 + temperature)))-0.6445)*(1-(373.16/(273.15 + temperature)))-1.976)*(1-(373.16/(273.15 + temperature)))+13.3185)*(1-(373.16/(273.15 + temperature)))))/1000

#wi,  leaf water vapor, converted to mole fraction
wi=ei/press_bar

#ea ambient vapor pressure (kPa)
ea_v <- (Relative_humidity/100)*ei

#wa, ambient water vapor, mole fraction
wa_v= ea_v/press_bar

# VPD in (kPa)
vpd =ei-ea_v 
round(vpd, digits=4 )

Here we can see the differences in Vapour presure defficity along these three times during the growing season.

#Total Conductance in mol
g = 1/(1/Stomatal_conductance + 1/br) # or 1/(rs+rb)

# E leaf transpiration (mol m-2 s-1) # barbour 2004
E_v_gs <- (wi-wa_v)*g

#ws, leaf surface water vapor, mole frac Roden 2000, Ball (1987),
ws_v_gs <- ((Stomatal_conductance*wi)-E_v_gs*(1-wi/2))/(Stomatal_conductance - E_v_gs/2)

# es vapor pressure at leaf surface (kPa) higher pressure  than the atmosphere
esl_v_gs <- ws_v_gs*press_bar

round(esl_v_gs, digits = 2)
round(ea_v, digits=2)

Lets recall some values that we already have estimated that will be important to keep track and estimate the expectation on the leaf water enrichment:

Rsource Source water
Rwater_vapor isotopic value of the water vapor
ei Water vapor presure saturated inside the leaf
ea_v Water vapor pressure atmosphere
esl_v_gs Water vapor pressure at leaf surface

And some constants like alpha_kappa, Kinetic Fractionation and alpha_kappa_beta, Kinetic Fractionation boundary layer

ak <- 1.0285    # alpha_kappa, Kinetic Fractionation 
akb <- 1.0189   # alpha_kappa_beta, Kinetic Fractionation boundary layer

# here we estimate the Leaf water enrichment:
Rleaf_v <- alpha_star*((ak*Rsource *
                              ((ei-esl_v_gs)/ei))+
                             (akb*Rsource *(esl_v_gs-ea_v)/ei)+
                             (Rwater_vapor*ea_v/ei))

# Rleaf_v now the leaf water is estimated in terms of R  and we must convert this value  on permill in rsmow Oxygen Isotope ratio from STANDARD SMOV
# Recall
rsmow <- 0.0020052
# we started on a value  of -10 permil now after  the evaporative enrichment  the leaf water has moverd  from -10 to 5 permil
leaf_water_18O <- ((Rleaf_v/rsmow)-1)*1000
round(leaf_water_18O, digits = 2)

PART 3: isotopic signal in the Cellulose

Isotopic fractionation during sugar synthesis and mixing during cellulose synthesis.

Now, after we know the water where the sugars are made and transported in the phloem, we can use the leaf water and the Source_18O to estimate the result on the cellulose, take into consideration that during cellulose synthesis, the source water will have some influence on the result, and we usually use a 42 percent exchange.

# Now   after we know the  water where the sugars are made and trasported in the floem  we can use the leaf water and the Source_18O to estimate the end result on the cellulose, take in consideration that during  cellulose sinthesisi  the source water will have some influence on the end result. and we usualy use a 42 percent exchange
# and some constants

epsilon_s <-  27 #Autotrophic fractionation Sucrose 
epsilon_c <- 27 # Autotrophic fractionation Cellulose ## 0.0084*Tmin.ts^2-0.51*Tmin.ts+ 33.172 
F_0 <- 0.42 # proportion of source water used during cellulose sinthesys

cellulose_d18O= F_0*(Source_18O + epsilon_c) + (1 - F_0) * (leaf_water_18O+ epsilon_s) 

PART 4: Comparison between modeled data and observations

Mechanistic modeling of tree‐ring cellulose δ18O and leaf water enrichment ΔL, and their sensitivities to relative humidity. (a) The Craig–Gordon model (CG) outputs are plotted from May to December, the period corresponding to cell wall thickening for all tree‐ring subdivisions. The yellow color lines correspond to the CG simulations of δ18O cellulose within the range of the proportion of exchangeable oxygen in cellulose (Pex) variability from 0.1 (thin line) to 0.7 (thick line). The gray horizontal bars correspond to the cell wall thickening (or cell‐maturation) period for each subdivision and are plotted on the mean of the δ18O measurements. The cell wall thickening phase was derived from microscopic phenological observations from Morino and Hughes (University of Arizona, unpublished observations, Table S3). The climate drivers used in the CG model to make δ18Ocel predictions were derived from PRISM data; (b) seasonal cycle of the CG mean leaf water enrichment (ΔL). The shaded gray color corresponds to the standard deviation of modelled ΔL. The dashed blue line in (a) and (b) indicates the average date of the monsoon onset; (c) modelled ΔL plotted versus relative humidity (RH) in % calculated at bimonthly resolution. The ΔL simulations were run using mean, minimum, and maximum RH. The range of RH for the hyper‐arid period of May and June are plotted in orange circles. The dark blue circles represent RH range observed for the monsoon season (July–August). The gray circles correspond to the range of RH observed from September to April (Belmecheri et al. 2018 PCE)

your estimations for the callulose oxygen Isotopes are:

round(cellulose_d18O, digits = 2)

your estimations for the leaf water oxygen Isotopes are:

round(leaf_water_18O, digits = 2)

Comparison between modeled data and observations,

Now with these values that you comared write down your observations differences and similarities between the modeled data in this excersise and the measurements shown in the figure 2 above.

How do you think this tool is useful?

What are the constrains or difficulties on implementing this model for reliable analysis?